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THERMAL  EXPANSION  OF  COMPOSITES  WITH  SPHERICAL 
AND  CYLINDRICAL  INCLUSIONS 


I.  INTRODUCTION 

A  composite  which  has  inclusions  of  different  thermal  expansivity  in  a 
matrix  will  develop  internal  stresses  as  the  temperature  is  changed.  The 
overall  coefficient  of  thermal  expansion  will  therefore  differ  from  the  simple 
volume  average. 

Previous  treatments  by  EshelbyO),  by  Wakashima  et  al.(2)  and  by 
Kerner(3)  used  the  stress  balance  at  the  interface  to  find  the  strain  field  in 
terms  of  the  misfit  between  inclusion  and  matrix,  and  are  confined  to  linear 
elasticity.  The  present  treatment  derives  these  relations  by  minimizing  the 
total  strain  energy.  It  uses  an  effective  medium  approximation,  i.e.  it 
replaced  the  properties  of  the  matrix  by  those  of  the  material  as  a  whole,  but 
assumes  that  the  short-range  shear  strain  about  each  inclusion  depends  only 
on  the  shear  modulus  of  the  matrix.  Also  it  expresses  the  result  in  terms  of 
bulk  moduli  and  shear  moduli,  since  this  is  convenient  when  separating  the 
strain  field  into  dilation  and  shear  components.  This  formulation  facilitates 
the  extension  of  the  theory  into  the  plastic  regime.  The  treatment  is 
confined  to  inclusions  which  are  spherical  or  randomly  oriented  long 
cyl inders. 
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II.  ELASTIC  FIELD  AROUND  A  SPHERICAL  INCLUSION 


Consider  a  spherical  hole  of  radius  R  cut  out  of  the  matrix,  with  a 
spherical  inclusion  of  radius  R  +  AR2  forced  into  the  hole  (if  AR2  >  0) ,  or 

placed  into  the  hole  with  the  surrounding  matrix  forced  to  join  the  inclusion 
(if  AR2  <  0).  The  hole  expands  or  contracts  to  a  radius  R  +  A  R. 

A  point  in  the  matrix,  of  distance  r  from  the  center,  undergoes  a 
displacement  u(r),  which  must  be  of  the  form 

u( r)  =  Ar  +  B/r2  ( 1) 

and  in  particular 

u(R)  =  AR  (2) 

The  spherical  inclusion  changes  its  radius  by  AR2  -AR;  it  changes  its 
volume,  and  there  is  an  associated  strain  energy.  If  the  fractional  volume  of 
inclusions  is  c,  the  energy  of  compression  of  the  inclusions,  per  unit  volume 
of  materi al  ,  i s 

Eincl  =  |  c  *1  (AR?  -  A  R)2/R2  (3) 

where  Kj  is  the  bulk  modulus  of  the  inclusion. 

The  term  Ar  in  (1)  corresponds  to  a  uniform  expansion  of  the  matrix;  the 
strain  energy  per  unit  volume  of  material  is 

^exp  =  j  9A2  (1-c)  (4) 

where  is  the  bulk  modulus  of  the  matrix. 

In  addition,  the  matrix  suffers  shear  because  of  the  non-di 1 atati onal 
displacement  u  =  B/r2.  The  principal  strains,  which  are  along  the  radial 
direction  and  any  two  mutually  perpendicular  directions  normal  to  it  are, 


'rr  =  3u/3r  =  -2B/r3 

(5a) 

yy  =  ezz  =  u/r  =  R/r3 

(5b) 

2 


^  v»  •/>  \r  ^  v*  ‘V  •-  - 


The  shear  strain  energy  density  becomes 

wsh  =  7  1  2^e i j 2  +  XA2 

=  viB2  r-«  (4+1+1)  =  6  yB2/r6  (6) 

since  the  dilation  A  =  err+eyy+ezz  vanishes.  Here  y  is  the  shear  modulus  of 
the  matrix,  \  is  the  other  Lame  modulus.  Integrating  over  the  volume  of  the 
matrix,  i.e.  from  r  =  R  to  infinity 

Eshear  =  7  4*r2  «sh(r)  <*r  =  (4ir/3)  6yB2/R3 
R 

=  (4wR3/3)  6y ( B/R3 ) 2  (7) 

Expressing  this  per  unit  volume  of  material,  where  c  can  be  equated  to  4wR3/3, 

Eshear  »  6rc  (8/R3)2  (8) 

The  total  strain  energy  per  unit  volume  of  material  thus  becomes 

E  =  (9/2)  c  Kj  (AR2  -  AR)2/R2  +  (9/2)(l-c)  KM  A2  +  fiyC  (B/R3)2  (9) 

where 

AR  =  R  (A  +  B/R3)  (10) 

Here  K j ,  Kfi  are  the  bulk  moduli  of  the  inclusion  and  the  matrix  respectively, 
u  is  the  shear  modulus  of  the  matrix,  c  is  the  volume  fraction  occupied  by 
inclusions,  A  and  B/R3  are  linear  strain  parameters  and  aR2  is  the  linear 
misfit  between  inclusion  and  matrix. 


One  can  define  two  strain  parameters 


0  =  B/R3 

(11) 

Y  =  AR2/R 

(1?) 

Now  the  condition  of  elastic  stability  is  that  the  total  strain  energy  E 


should  be  a  minimum;  two  stability  conditions  must  be  satisfied 


3E/9A  =  0  (13a. 

and 

3E/30  =  0  ,  (13b) 


Since  (a(?2  -  AR)/R  =  y-A-0,  equation  (13a)  yields 
9c  Kj  (y-A-S)  *  9( 1-c)  KMA 
and  ( 13b)  yields 


9c  K j  (y-A-6)  =  12  uc6 

Eliminating  the  left-hand  side  of  both  equations 


c  B  4u_ 
1-c  6  IRm 


(14a) 

(14b) 


(IB) 


Eliminating  A  from  (14a  and  b) 


8 


Y 


4uc  i-l 
3KM(T-c)  J 


and  eliminating  0  from  (15)  and  (16) 


A 


c  4u 
Y  1-c  3Km 


[1 


4u  4pc  i-l 
3K  j  3KM( 1-C )  1 


(16) 


(17) 


III.  APPLICATION  TO  THERMAL  EXPANSION 

The  volume  expansion  comes  from  two  sources:  the  A-field  is  a  uniform 
expansion  of  the  matrix  and  the  included  cavities,  while  the  0-field 
represents  an  additional  expansion  of  the  included  cavities.  This  latter 
field  is  non-di 1 ati onal  ,  and  transmits  that  expansion  to  the  outer  boundaries, 
since  =  4itr2B/r2  is  independent  of  r.  Therefore 


6 V/V  =  3A  +  3c6 


(1ft) 
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Now  let  the  matrix  have  a  volumetric  thermal  expansion  coefficient 
dAM/dT,  and  the  inclusions  a  coefficient  dAj/dT,  then  the  overall  coefficient 
of  thermal  expansion  is 

dA/dT  =  dAM/dT  +  ^  ( 6V/V) 

=  dAM/dT  +  (3A  +  3c6)  (19) 

where  A  and  g  can  be  expressed  in  terms  of  y  through  (16)  and  (17). 

Now  for  an  isolated  inclusion 

3  dy/dT  =  dAi/dT  -  dAM/dT  (20) 

but  if  c  is  not  small,  dAM/dT  should  be  replaced  by  the  actual  expansion  of 
the  composite,  which  is  dA/dT.  This  takes  account,  to  some  degree,  the 
interaction  between  inclusions;  in  addition.  Km  should  be  replaced  by  the 
average  bulk  modulus  of  the  material 


KM  =  (1-c)  Km  +  c 


On  the  other  hand,  the  shear  field  B/r3  is  short-range,  and  if  the  inclusions 
do  not  touch  too  frequently,  u  should  be  taken  as  the  shear  modulus  of  the 
matrix.  Thus  replacing  (20)  by 

3dy/dT  =  dA  j/dT  -  dA/dT  (22) 

equation  (19)  becomes,  using  (16)  and  (17) 

dA/dT  =  dAM/dT  +  3c  (dy/dT)  F(u)  (23) 


Using  (22),  one  finally  obtains 


where 


( 1+cF )  dA/dT  =  dA m/ d T  +  cF  dAj/dT 


1  +  4m/3(1-c)Km 
1  +  4u/3Kj  +  4cy/3( 1-c)Km 


is  a  function  of  u. 


To  obtain  the  net  dilation  over  a  finite  temperature  interval,  one  has  to 
integrate  (24)  over  temperature,  i.e. 

A(T)  «  /  dT  (26) 

To 

It  is  relatively  straight  forward  to  take  account  of  temperature  dependences 
of  dAfi/dT  and  dAj/dT.  Similarly  Kj  and  can  be  regarded  as  temperature 
dependent,  although  usually  only  weakly  so. 

The  shear  modulus  u  is  defined  by 


2p  =  do-j/dej  i*j  (27) 

where  the  stress  and  strain  tensors  a  and  e  are  expressed  in  the  contracted 
notation.  Now  u  is  not  only  a  function  of  temperature,  but  in  the  plastic 
regime  also  a  function  of  prior  shear  strain.  Since  the  shear  strain  is  a 
function  of  radial  distance  r  (see  equations  5  and  6),  u  is  also  a  function  of 
r,  and  therefore  the  present  treatment  is  strictly  speaking  invalid  once  the 
plastic  regime  is  reached. 

IV.  FFF1CT  OF  PLASTICITY 

Although  p,  which  depends  on  prior  shear  strain,  is  different  at 
different  points  of  the  matrix,  since  the  shear  strain  is  a  function  of 
position,  one  notes  that  most  of  the  shear  strain  energy  of  (7)  resides  in  the 
immediate  vicinity  of  the  inclusion,  where  the  shear  strain  is  a  maximum,  and 
is  given  in  magnitude  by  S  =  B/R3 .  To  a  good  approximation  one  may  thus  regard 
u  of  equation  (24)  to  be  only  a  function  of  6  and  T. 
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model  of  plasticity.  As  T  is  increased  by  a  finite  step  dT,  the  incremental 
expansion  dA  must  be  determined  from  (24),  and  the  increment  in  shear  strain 
then  determined  from 

d0  =  (1/3)  (dA!  -  dA  )  /  f  ( y )  (31) 

The  new  value  of  A,  namely  A+dA,  and  the  new  value  of  3,  namely  3+d3,  are  then 
used  to  determine  a  new  value  of  y,  and  this  value  is  used  in  (34)  and  (31) 
for  the  next  temperature  increment.  This  procedure  seems  suited  for  a 
programmable  computer, 

V.  ELASTIC  FIELD  AROUND  AN  INFINITE  CYLINDER 

The  displacement  field  u(r,o,z)  about  an  infinite  cylinder  obeys  the 
elastic  equations 

jfijf  +  7)  =  0,  a2u/az2  =  0  (32) 

and  the  solution  is  of  the  form 

ur  =  ( Ar  +  B/r)  ,  uz  =  e0z  (33) 
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The  principal  strains  are 
ezz  =  eo 

err  =  3u/3r  =  A  -  B/r2 
e00  =  u/r  =  A  +  B/r2 


and  the  dilation  is 


A  =  2A  +  er 


The  principal  stresses  are 


(34a) 

(34b) 

(34c) 


arr  =  2y(A  -  B/r2)  +  AA 
O0Q  =  2y(A  +  B/r2  +  AA 

Ozz  =  2yeQ  AA 


(36a) 

(36b) 

(36c) 


where  A,  y  are  the  Lame  moduli.  The  strain  energy  density  is  given  by 

W(r)  =  1/2  orrerr  +  1/2  o00  ege  +  1/2  ozzezz 

=  y( A-B/r2)2  +  y(A+B/r2)2  +  ue02  +  1/2  AA2 

=  2yA2  +  2ul32/r4  +  Me02  +  1/2  A(2A+e0)2  (37) 

Now  consider  a  random  assembly  of  cylinders,  filling  a  fraction  c  of  the 
entire  volume.  Let  the  length  of  each  cylinder  experience  the  same  fractional 
increase  as  any  random  line  in  the  matrix,  so  that  e0=A.  The  strain  energy  of 
the  matrix  is  obtained  by  integrating  the  B-field  energy  density  from  r=R  to 
r=®,  and  adding  to  it  the  dilatational  energy  of  the  A-field,  which  is  a 
uniform  expansion.  Thus 


Eb  =  2uB2  /  2*rr-4  dr  =  2ttuB2/R2 
R 

=  2u(  B2/R4  )  (ttR2) 


Note  that  nR2  is  the  volume  per  unit  length  of  the  inclusions,  so  that,  per 
unit  volume,  Eb=2|/(B2/R4  )c .  For  the  shear  modulus  of  the  matrix  we  take  that 
of  the  matrix  material,  assuming  that  since  most  of  the  strain  energy  is  in 
the  immediate  vicinity  of  the  inclusion,  the  significant  material  is  that  of 
the  matrix.  Now  K^= ( A+2y/3)  is  the  hulk  modulus  of  the  matrix.  The 
dilational  energy  can  thus  also  be  written  as  1/2  Km(3A)2,  and  since  the 
matrix  occupies  a  fractional  volume  (1-c)  only,  this  energy  must  be  multiplied 
by  that  factor  (1-c).  Furthermore,  should  really  be  the  volume  average  of 
the  matrix  and  the  inclusions,  so  that,  per  unit  volume 


EA  =  1/2  ( 1-c) ( 3A)2  [(1-c)  KM  +  c  K 1} 


where  Kj  is  the  bulk  modulus  of  the  inclusion  material.  The  strain  energy  of 
the  matrix  material,  per  unit  volume  becomes  thus 


Em  =  EA  +  EB  =  1/2  ( 1-c) Kw  ( 3A) 2  +  2cpMB2/R4 


where  is  the  volume-averaged  bulk  modulus. 

Now  the  strain  field  arises,  in  the  first  place,  because  each  cylindrical 
inclusion,  when  unstrained,  occupies  a  different  volume  than  the  cylindrical 
cavity  of  the  unstrained  matrix.  Let  the  original  radius  of  the  cyliodircal 
inclusion  be  R+aR?;  it  sits  in  a  matrix  cavity  of  radius  R+AR,  where 


AR  =  AR  +  B/R 


Let  us  also  assume  that  originally  there  is  a  longitudinal  strain  between  the 
unstrained  matrix  and  the  inclusion,  of  magnitude 


e 2  -  AR2/R 


I  SI  *•  I  " if 


KSS 


w? 


th  ,.i  rr 


Thus  the  original  length  of  the  inclusion 

l  +  e2  =  l  +  AR2/R  (42" ) 

is  changed  to  1+A.  The  inclusion  thus  has  the  following  principal  strains: 

AR2/R  -  aR/R  (2  transverse  components) 

AR2/R  -  A  (1  longitudinal  component) 

and  a  net  dilation 

A  i  =  3A R2/R  -  2AR/R  -  A  (43a) 

and  in  view  of  (41) 

Aj  =  3A R2/R  -  3A  -  2B/R2  (43b) 

This  strain  field,  uniform  within  the  inclusion,  can  be  resolved  into  a 
dilation  a  of  (43b),  and  three  principal  strains  of  zero  dilation,  of 
magni tude 

-2B/R2,  B/R2  and  B/R2 

so  that,  per  unit  volume  of  material,  the  strain  energy  within  the  inclusions 
becomes 

Ei  =  1/2  cK[A 2  +  6cp jB2/R4  (44) 

while  the  total  strain  energy  per  unit  volume  becomes 

E  =  EM  +  =  (9/2)  Km  ( 1-c)  A2  +  ?cuf1  B2  + 

+  6cu i  32  +  ( 9/ 2 ) c  Kr  (Y-A-20/3)2  (45) 


where  3=R/R2  and  y=AR2/R. 

The  conditions  of  stability  are  that 

9  E/3  A  =  0  and  3  E/ 3  0  =  0 

which  becomes 

9 Km  (1-c) A  =  9Kj  c  ( y-A-20/3) 
4(um+3ui)cB  =  6K j  c  (y-A-2S/3) 


(47a) 

(47b) 


Tt 

1 

I 

lV.i 

l;?; 


IS 


Eliminating  y  we  obtain 

A  =  -|  -y^r  (um+3hi)  6/Km 

Substituting  (48)  into  (47b)  we  obtain 

3=2  Y  [i  +  EMl2UL  +  _L_  iiMilHl]' 
2  Kj  1-c  Km 

Finally,  from  (48)  and  (49) 


A  =  Y  [WM+frl]  [j  +  u_m+3ui  +  _c_  um+3ui  j* 
1-C  Km  Kj  1-C  Km 

+  ♦  c/d-or1 


Thus,  given  y,  one  can  find  A  and  g. 


VI.  VOLUME  CHANGE  AND  THERMAL  EXPANSION 


The  strain  field  just  discussed  causes  a  net  dilation 

A  =  3A  +  2cB  (51) 

The  first  term  is  a  uniform  expansion  of  the  matrix  and  of  the  cylindrical 
cavities  which  hold  the  inclusions;  the  second  term  is  the  additional 
expansion  which  does  not  dilate  the  matrix  but  transmits  the  expansion  of  the 
cylinders  to  the  outside  surface.  This  can  be  expressed  in  terns  of  the 
misfit  parameters  y  through  (49)  and  (50).  Writing 

um+3ut  c  pm+3uj 

■  1  *  m  *  l-c  k;  (62) 

equation  (51)  becomes 

pm+3ut 

4  “  3c*  O  +  (NcW  f'“  <53> 

In  the  special  case  when  Kj=Km  and  hence  Kj=Km  ,  (53)  becomes  simply 
A=3cy,  so  that  the  excess  volume  of  the  cylindrical  inclusions,  weighted  by  c, 
becomes  simply  the  net  increase  in  volume:  the  compression  of  the  inclusions 


..frrtV.n 


by  the  matrix  is  compensated  by  the  expansion  of  the  matrix.  In  general, 
however,  the  expansion  is  larger  than  3cy  if  Kj  >  Km,  and  vice  versa. 

Suppose  now  that  the  misfit  parameter  y  is  caused  by  a  difference  in 
thermal  expansivity  of  matrix  and  inclusions,  so  that  as  in  (20) 

3dy/dT  =  dA  j/dT  -  dAM/dT  (54) 

Writi ng 

F(uM»mi)  =  (1  +  -UJi+~Uf)/f  (55) 

(l-c)KM 

where  f  is  given  by  (52),  adding  the  thermal  expansion  of  the  matrix  to  that 
due  to  the  misfit  (equation  53)  we  obtain 

dA/dT  =  dAffiT  +  3cF  dy/dT  (56) 

With  the  same  self-consistency  approximation  as  was  made  for  spherical 
inclusions  we  replace  dAM/dT  in  (54)  -but  not  in  the  first  term  of  (56)  -by 
dA/dT,  and  obtain,  as  in  (24),  that 

(1+cF)  dA/dT  =  dAM/dT  +  cF  dAi/dT  (57) 

the  only  difference  being  that  F,  previously  defined  by  (25),  is  now  given  by 
(55). 

When  accounting  for  plasticity,  we  must  remember  that  F  is  now  a  function 
of  two  shear  moduli,  pm  an(1  mi»  We  must  still  make  the  approximation,  not  as 
well  justified  for  cylinders  as  for  spheres,  that  pm  is  just  a  function  of  6 
and  T,  since  the  shear  strain  energy  in  the  matrix  resides  mainly  just  near 
the  interface.  In  the  inclusions,  however,  the  shear  strain  in  uniform,  so 
that  no  approximation  is  incurred  when  treating  pj  as  just  a  function  of  p  and 
T.  One  can  then  proceed  analogously  to  the  case  of  spherical  inclusions, 
obtaining  A (T)  by  numerical  integration. 
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VII.  CONCLUSION 

By  minimizing  the  strain  energy  and  making  use  of  the  fact  that  the 
strain  field  is  composed  of  a  uniform  dilation  of  the  matrix  and  a  short-range 
shear  strain  field,  it  is  possible  to  obtain  closed  expressions  for  the 
overall  thermal  expansion  in  the  case  of  spherical  inclusions  and  of  randomly 
oriented  long  cylindrical  inclusions.  A  procedure  was  obtained  for  extending 
the  results  to  the  case  of  plastic  flow  of  the  matrix  about  the  inclusion,  if 
plastic  behavior  can  be  expressed  as  an  effective  shear  modulus  as  function  of 
temperature  and  shear  strain. 

Expressions  were  also  obtained  for  the  dilation  of  the  matrix  in  terms  of 
the  misfit  between  matrix  and  inclusion,  and  thus  also  in  terms  of  the  overall 
expansion.  These  could  be  used  to  test  the  model  by  comparing  dilatometry 
with  changes  of  the  lattice  spacing  of  the  matrix.  These  results  could  also 
be  used  in  related  problems,  such  as  comparing  the  swelling  due  to  radiation 


damage  to  the  corresponding  changes  of  the  lattice  spacing  of  a  material. 
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